*! Inspirit of Lian Yujun professor of Sun Yat-Sen University.

*! Authors:       Shi Jinyu(18333158)
*! About Authors: Undergraduates from Department of Finance, Lingnan College, 
*!                Sun Yat-Sen University.
*! Reference:     Thanks greatly for the code offered by Professor Yujun Lian(Arlion), 
*!                Department of Finance, Lingnan College, Sun Yat-Sen University.
*! Date:          2021.1.31	



	global path "D:\stata15\ado\personal\FE_Bech\EfficientFrontier"
	cd $path

	cap mkdir data
    cap mkdir refs
    cap mkdir out
    cap mkdir adofiles  
	
	global D    "$path\data"      //范例数据
	global R    "$path\refs"      //参考文献
	global out  "$path\out"       //图形和表格结果
	adopath +   "$path\adofiles"  //添加外部命令所在路径  
	cd $D
	set scheme s2color //设置彩色绘图模板
  

*------------------
*- 0 安装外部命令
*------------------

  *-设定安装外部路径的位置
    net set ado "$path\adofiles"
    net set other "$D"  //外部命令附带的 example datasets 等存放于此
	
  *-安装外部命令
    ssc install cntrade, replace //下载个股和市场指数实时数据, SJ 14(2):381--388
	ssc install openall, replace //批量合并
	ssc install mvport, replace //投资组合模型的估计程序


*---------------------
*- 2.1 下载个股交易资料
*---------------------
  *-下载并保存2020年整年科大讯飞(002230)、贵州茅台(600519)和药明康德(603259)///
  *-的股票交易数据
  
  local list "002230 600519 603259 "
  cntrade `list'
  openall `list'
  gen year = year(date)
  keep if year == 2020
  save "mstocks_long.dta", replace //保存合并后的数据
  
  
 *-纵横变换: reshape 
    use "mstocks_long.dta", clear          //long format
    keep stkcd date clsprc //仅保留需要纵横变换的变量
    xtset date stkcd         
    reshape wide clsprc, i(date) j(stkcd)  
  save "mstocks_wide.dta", replace  //wide format
  
  use "mstocks_long.dta", clear          //long format
    keep stkcd date clsprc //仅保留需要纵横变换的变量
    xtset date stkcd           //这一步很重要
    reshape wide clsprc, i(date) j(stkcd)  
  save "mstocks_wide.dta", replace  //wide format
*---------------------
*- 2.2 观察股票走势
*---------------------
  /* 股价标准化*/
  use mstocks_wide, clear
  foreach i of varlist clsprc*{
    if "`i'" != "date"{
        replace `i' = (`i' / `=`i'[1]') * 100
        format `i' %6.2f
      }
   }

  tw
  line  clsprc*   date, ///
    xla(21915(45)22281, ang(20)) ///
    leg(order(1 "科大讯飞" 2 "贵州茅台" ///
        3 "药明康德" ) position(1) ring(0)) ///
	title("股价走势图（2020年）") ///
	xtitle("日期") ///
	ytitle("标准化股价") 
  graph export "$out\price_trend.png", replace   
  
*---------------------
*- 2.3 计算收益率与协方差
*---------------------	  
  use mstocks_wide, clear
   //由于股票交易日之间可能存在间隔，所以不能直接用日期作为时间变量
  gen dateid = _n
  tsset dateid
  foreach i of varlist clsprc*{
        gen l`i' = l.`i'
        replace `i' = log(`i'/l`i')
        drop l`i'
  } 
  drop dateid
  rename (clsprc2230 clsprc600519 clsprc603259) (科大讯飞 贵州茅台 药明康德)
  list in 1/5	
  save return, replace
  
  use return, clear
  drop in 1
  local j = 1
  foreach i of varlist 科大讯飞 贵州茅台 药明康德{
        qui sum `i'
        local r`j' = r(mean)
        local j = `j' + 1
  }
  di "`r1', `r2', `r3'"	
  mat rets = (`r1', `r2', `r3')	
	
  corr 科大讯飞 贵州茅台 药明康德, cov
  mat cov = r(C)
*---------------------
*- 3 构造资产组合
*---------------------	 
  clear 
  set obs 11
  qui egen double w1=fill(0(0.1)1)
  qui format w1 %tg
  qui gen w2=1-w1
  qui format w2 %tg

  qui gen r_P = w1*rets[1,1] + w2*rets[1,2]
  qui gen varP= w1^2*cov[1,1] + w2^2*cov[2,2] + 2*w1*w2*cov[2,1]
  qui gen sdP=sqrt(varP)
  twoway (scatter r_P sdP, msize(tiny) mlabel(w1) mlabcolor(edkblue)), ///
  ytitle(收益率) ylabel(#5) xtitle(波动性) /// 
  title(Frontier using 2 stocks(科大讯飞 & 贵州茅台))
  graph export "$out\Frontier_using_2_stocks.png", as(png) replace



  clear 
  set obs 101
  qui egen double w1=fill(0(0.01)1)
  qui format w1 %tg
  qui egen double w2=fill(0(0.01)1)
  qui gen w3=1-w1-w2
  qui format w2 %tg
  qui format w3 %tg

  qui gen r_P = w1*rets[1,1] + w2*rets[1,2]+w3*rets[1,3]
  qui gen varP= w1^2*cov[1,1] + w2^2*cov[2,2] +w3^2*cov[3,3]+2*w1*w3*cov[1,3] ///
   +2*w2*w3*cov[2,3]+2*w1*w2*cov[2,1]
  qui gen sdP=sqrt(varP)
  twoway (scatter r_P sdP, msize(tiny) mlabcolor(edkblue)), ///
    ytitle(收益率) ylabel(#5) xtitle(波动性) ///
    title(Frontier using 3 stocks(科大讯飞 & 贵州茅台 & 药明康德))
  graph export "$out\Frontier_using_3_stocks.png", as(png) replace

*---------------------
*- 4.1 蒙特卡洛模拟
*---------------------	
cap prog drop front
prog def front, rclass
    version 15.0
    use return, clear
    local w1 = runiform()
    local w2 = runiform()
    local w3 = runiform()
    mat weight = (`w1' \ `w2' \ `w3' )
    mat list weight
    mat weight = weight / (`w1'+`w2'+`w3')
    mat list weight
    ret scalar w1 = weight[1, 1]
    ret scalar w2 = weight[2, 1]
    ret scalar w3 = weight[3, 1]
    local j = 1
    foreach i of varlist _all{
        if "`i'" != "date"{
            qui sum `i'
            local r`j' = r(mean)
            local j = `j' + 1
        }
    }
    mat rets = (`r1', `r2', `r3')
    mat a =  rets * weight
    mat list a
    ret scalar ret = a[1, 1]
    corr 科大讯飞 贵州茅台 药明康德, cov
    ret list
    mat cov = r(C)
    mat b =  weight' * cov * weight
    mat list b
    ret scalar var = b[1, 1]
    ret scalar std = sqrt(b[1, 1])
end

simulate ret = r(ret) var = r(var) std = r(std) w1 = r(w1) w2 = r(w2) ///
  w3 = r(w3) , reps(5000): front
  save dataset, replace
tw ///
sc ret std, msize(*0.01) xtitle(波动率) ytitle(收益率) ///
    title("Portfolio") msymbol(o) ///
    xla(#6, format(%6.3f)) yla(, format(%6.4f))
graph export "$out\Portfolio.png", as(png) replace
*---------------------
*- 4.2 计算并画出有效边界
*---------------------	
mata:
 mata clear
  void min_std(real scalar todo, real vector w12, std, g, H)
  {
    real scalar w3
    w3 = 1 - sum(w12)
    real vector r_mean
	r_mean= st_matrix("rets") //将储存于 Stata 的矩阵导入 Mata环境
    real vector w
    w = (w12, w3)
	
	real scalar r_P
    r_P = w * r_mean'
    real matrix var_P
    var_P = st_matrix("cov")
    var_P = w * var_P * w'
    std = -sqrt(var_P[1, 1])
  }

    S = optimize_init()
    optimize_init_evaluator(S, &min_std())
    optimize_init_params(S, J(1, 2, 0))
    wh = optimize(S)

    w = (wh, 1-sum(wh))
    r_Pm = (wh, 1-sum(wh)) * st_matrix("rets")'
    var_P = (w * st_matrix("cov") * w')
    std = sqrt(var_P[1, 1])

   w

   r_Pm
 
   std

end

mata:
 mata clear
  void min_std(real scalar todo, real vector w123, std, g, H)
  {
    real vector rmean
    rmean = st_matrix("rets")
    real scalar rp
    real vector w
    w = w123
    rp = w * rmean'
    real matrix variancep
    variancep = st_matrix("cov")
    variancep = w * variancep * w'
    std = -sqrt(variancep[1, 1])
  }
  
  std_vector = J(1, 40, 0)
  j=1
//最小化方差对应的收益率为0.00196，最大收益率为仅购入贵州茅台时的0.00235
    for (i = 0.00196; i < 0.00235; i = i + 0.00001) 
    {
   S = optimize_init()
      optimize_init_evaluator(S, &min_std())
      optimize_init_technique(S, "nr")
      optimize_init_constraints(S, ((st_matrix("rets")\ J(1, 3, 1)), (i \ 1)))
      optimize_init_params(S, J(1, 3, 0))
    wh = optimize(S)

    variancep = (wh * st_matrix("cov") * wh')
    std = sqrt(variancep[1, 1])
    std_vector[1, j] = std
    j = j + 1
   }

  st_matrix("std_vector", std_vector)

end


clear
set obs 40
gen ret = (_n + 195) / 100000
gen std = .
forval i = 1/40{
    replace std = std_vector[1, `i'] in `i'
}
gen front = 1
drop if std == 0
save front, replace


tw ///
sc ret std if front == 1,  xtitle(波动率) ytitle(收益率) ///
    title("Efficient Frontier using 3 stocks") ///
	 msize(tiny) msymbol(D) mc("blue") leg(off) || ///
	scatteri .0019607171   .0171638423 "最小方差点", ///
    mc("red")  msymbol(o) mlabc("pink")  ///
 
graph export "$out\Efficient_Frontier_using_3_stocks.png", as(png) replace

use dataset, clear
append using front
tw ///
sc ret std, msize(*0.01) xtitle(波动率) ytitle(收益率) ///
    title("Efficient Frontier using 3 stocks") msymbol(o) ///
    xla(#6, format(%6.3f)) yla(, format(%6.4f))|| ///
scatteri  .0019607171   .01716384 "最小方差点", ///
    mc("pink")  msymbol(o) mlabc("pink") || ///
sc ret std if front == 1, msize(tiny) msymbol(D) ///
    mc("blue") leg(off)
graph export "$out\Investment_portfolio_using_3_stocks.png", as(png) replace
 
*---------------------
*- 5 快捷命令
*---------------------	

use return.dta, clear
gmvport 科大讯飞 贵州茅台 药明康德, noshort
 
use return.dta, clear
qui efrontier 科大讯飞 贵州茅台 药明康德
graph export "$out\EFrontier.png", as(png) replace
 
